###################################################################
### Competing principals model, 20 chains on 10 samples, Mexico ###
###################################################################

rm (list=ls())
print (version)

# platform       x86_64-pc-linux-gnu         
# arch           x86_64                      
# os             linux-gnu                   
# system         x86_64, linux-gnu           
# status                                     
# major          4                           
# minor          3.0                         
# year           2023                        
# month          04                          
# day            21                          
# svn rev        84292                       
# language       R                           
# version.string R version 4.3.0 (2023-04-21)
# nickname       Already Tomorrow            

library (runjags)
library (coda)
library (reshape2)
library (random)
library (scales)

# Set working directory
savePathRep <- "/all_data/"
setwd (savePathRep)


all.seeds <- matrix (data= c(16353, 1139, 4301, 6615,
                             4843, 1608, 1617, 1822,
                             1939, 3549, 86687, 1267,
                             3035, 26666, 6630, 8199,
                             8261, 5051, 8882006, 5213),
                     ncol=2, byrow = FALSE)

row <- c(1:10)

for (j in 1:2){
for (i in 1:10){
  which.sample <- row[i]  # for sample_mean, use "mean"
  which.process <- ifelse (j==1, "a", "b")
  this.seed <- all.seeds[which.sample, j]
  # for sample_mean, use this.seed <- 1972 for a or 1973 for b
  
  load (paste0 ("mexico_rollcall_long_format_", which.sample, ".RData")) 

  load ("mexico_initial_values.Rda")

  potential.fix.bill <- c(22,89,144,150,174,217)

  pressure <- function(x1, x2, x3) { 
    k <- rep(0, length(x1))
    for (i in 1:length(x1)){
      k[i] <- ifelse(max(x1[i], x2[i], x3[i])>0, max(x1[i], x2[i], x3[i]), min(x1[i], x2[i], x3[i]))
    }
    k
  }


  compete.bis.v="model { 
     for (i in 1:n.obs) {
       probit(CJR[i]) <- beta[vote[i]]*theta[dep[i]] + eta[vote[i]]*delta[dep[i]] - alpha[vote[i]] + 
                        lambda[1]*P1[obs[i]] + lambda[2]*P2[obs[i]] + lambda[3]*P3[obs[i]] + 
                        lambda[4]*P4[obs[i]] + lambda[5]*P5[obs[i]] + lambda[6]*P6[obs[i]] + 
                        lambda[7]*P7[obs[i]] + lambda[8]*P8[obs[i]] + lambda[9]*P9[obs[i]] +
                        lambda[10]*P10[obs[i]] + lambda[11]*P11[obs[i]] + lambda[12]*P12[obs[i]] +
                        lambda[13]*P13[obs[i]] + lambda[14]*P14[obs[i]] + lambda[15]*P15[obs[i]] +
                        lambda[16]*P16[obs[i]] + lambda[17]*P17[obs[i]] + lambda[18]*P18[obs[i]] +
                        lambda[19]*P19[obs[i]] + lambda[20]*P20[obs[i]] + lambda[21]*P21[obs[i]]
         y[i] ~ dcat(Q[i,1:3]);         # code Nay as 1, NA as 2, Aye as 3
         Q[i,1] <- pow(probVoteDisagree*(1-CJR[i]), p1[i]) * pow(probVoteAgree*(1-CJR[i]), p2[i]) * pow(probVoteUncertain*(1-CJR[i]), (1-p1[i]-p2[i]))
         Q[i,2] <- 1 - Q[i,1] - Q[i,3]
         Q[i,3] <- pow(probVoteAgree*CJR[i], p1[i]) * pow(probVoteDisagree*CJR[i], p2[i]) * pow(probVoteUncertain*CJR[i], (1-p1[i]-p2[i]))
   }
   # PRIORS
   probVoteDisagree ~ dbeta(1,1)  # This is a flat prior
   probVoteAgree ~ dbeta(1,1)     # This is a flat prior
   probVoteUncertain ~ dbeta(1,1) # This is a flat prior
   for(j in 1:21) { lambda[j] ~ dnorm(0,1) } 
   for(j in 1:n.item) { alpha[j] ~ dnorm(0,1) }
  # Alternative identification: Beta (discrimination, dimension 1)
  for(j in 1:(fix.bill[1]-1)) { beta[j]  ~ dnorm(0,1) }
  beta[fix.bill[1]]~ dnorm(0,1)
  for(j in (fix.bill[1]+1):(fix.bill[2]-1)) { beta[j]  ~ dnorm(0,1) }
  beta[fix.bill[2]] <- 1.5
  for(j in (fix.bill[2]+1):(fix.bill[3]-1)) { beta[j]  ~ dnorm(0,1) }
  beta[fix.bill[3]] ~ dnorm(0,1)  #<- 2
  for(j in (fix.bill[3]+1):(fix.bill[4]-1)) { beta[j]  ~ dnorm(0,1) }
  beta[fix.bill[4]] <- -1
  for(j in (fix.bill[4]+1):(fix.bill[5]-1)) { beta[j]  ~ dnorm(0,1) }
  beta[fix.bill[5]] <-  2.5
  for(j in (fix.bill[5]+1):(fix.bill[6]-1)) { beta[j]  ~ dnorm(0,1) }
  beta[fix.bill[6]] ~ dnorm(0,1)
  for(j in (fix.bill[6]+1):n.item) { beta[j]  ~ dnorm(0,1) }
  # Eta (discrimination, dimension 2)
  for(j in 1:(fix.bill[1]-1)) { eta[j]  ~ dnorm(0,1) }
  eta[fix.bill[1]] ~ dnorm(0,1)
  for(j in (fix.bill[1]+1):(fix.bill[2]-1)) { eta[j]  ~ dnorm(0,1) }
  eta[fix.bill[2]] <- 2
  for(j in (fix.bill[2]+1):(fix.bill[3]-1)) { eta[j]  ~ dnorm(0,1) }
  eta[fix.bill[3]] ~ dnorm(0,1) #
  for(j in (fix.bill[3]+1):(fix.bill[4]-1)) { eta[j]  ~ dnorm(0,1) }
  eta[fix.bill[4]] <- -2 #
  for(j in (fix.bill[4]+1):(fix.bill[5]-1)) { eta[j]  ~ dnorm(0,1) }
  eta[fix.bill[5]] <- 0 #
  for(j in (fix.bill[5]+1):(fix.bill[6]-1)) { eta[j]  ~ dnorm(0,1) }
  eta[fix.bill[6]] ~ dnorm(0,1) #
  for(j in (fix.bill[6]+1):n.item) { eta[j]  ~ dnorm(0,1) }

  # ideal points
  for(i in 1:n.legs)  { theta[i] ~ dnorm(0,1) }
  for(i in 1:n.legs)  { delta[i] ~ dnorm(0,1) }
  ## Notes
    # n.party.legislators, n.legislators, and n.item are data (counters for the size of vectors)
    # y, party, party.position, pressure, and R are all data:
    # y is an i*j rollcall matrix with 1 (Aye), 0 (Nay), and NA (unobserved)
    # party is a vector of length i with 1=PRI, 2=PAN, 3=PRD
    # party.position is a matrix of size 3*j with the parties' official positions (1 or 0; NA not allowed at this point)
    # pressure is a matrix of size 3*j with continuous scores for each party's pressure on vote j
    # R is an i*j indicator matrix for legislators/votes for which pressure is expected to be effective
  }"

position.vec<- position.vec.strict

compete.data <- dump.format(list(y=car::recode (good.RC$RC, "1=3; 0=1; NA=2")
                                 , p1=ifelse(position.vec==1, 1, 0)
                                 , p2=ifelse(position.vec== -1, 1, 0)
                                 , n.obs=nrow(good.RC)
                                 , obs=seq(1:nrow(good.RC))
                                 , vote=as.numeric(good.comp.party$vote)
                                 , dep=good.comp.party$final.leg.no
                                 , n.legs=max(good.comp.party$final.leg.no)
                                 , n.item=max(as.numeric(good.comp.party$vote))
                                 , P1=good.R7$pressure #PRI
                                 , P2=good.R8$pressure #PAN
                                 , P3=good.R9$pressure #PRD
                                 , P4=pressure(good.R1$pressure, good.R2$pressure, good.R3$pressure)
                                 ,P5=good.R67$pressure   #PR
                                 ,P6=pressure(good.R27$pressure, good.R26$pressure, good.R28$pressure)
                                 ,P7=good.R64$pressure  #poverty dummy
                                 , P8=pressure(good.R69$pressure, good.R68$pressure, good.R70$pressure)
                                 ,P9=good.R65$pressure    # year
                                 ,P10=pressure(good.R30$pressure,good.R29$pressure, good.R31$pressure)
                                 ,P11=good.R75$pressure   #victory margin dummy
                                 ,P12=pressure(good.R51$pressure, good.R50$pressure,good.R52$pressure)
                                 ,P13=good.R72$pressure  #sta.muni.leg
                                 ,P14=pressure(good.R48$pressure, good.R47$pressure, good.R49$pressure)
                                 ,P15=good.R77$pressure  #CabinetExp
                                 ,P16=pressure(good.R37$pressure, good.R33$pressure, good.R41$pressure)
                                 ,P17=good.R71$pressure   #prezposition
                                 ,P18=pressure(good.R45$pressure, good.R44$pressure, good.R46$pressure)
                                 ,P19=good.R73$pressure    #commPres
                                 ,P20=pressure(good.R36$pressure, good.R32$pressure, good.R40$pressure)
                                 ,P21=good.R56$pressure    #pcbrequest
                                 , fix.bill=potential.fix.bill))

compete.parameters = c("theta", "delta", "alpha", "beta", "eta", "lambda", "deviance"
                       , "probVoteDisagree", "probVoteAgree", "probVoteUncertain")


##Initials for alternative identification
eta1 = initials$eta + runif(length(initials$eta), -0.2, 0.2)
beta1 = initials$beta + runif(length(initials$beta), -0.2, 0.2)
eta1[potential.fix.bill[4:5]] <- c(NA,NA)
eta1[potential.fix.bill[2]] <- c(NA)
beta1[potential.fix.bill[2]] <- c(NA)
beta1[potential.fix.bill[4:5]] <- c(NA,NA)

compete.inits.1 <- function (){
  dump.format (
    list(
      theta = initials$theta + runif(length(initials$theta), -0.2, 0.2),
      delta = initials$delta + runif(length(initials$delta), -0.2, 0.2),
      alpha = initials$alpha + runif(length(initials$alpha), -0.2, 0.2),
      eta = eta1,
      beta = beta1
         , probVoteAgree = rbeta (1,1,1)
         , probVoteDisagree = rbeta (1,1,1)
         , probVoteUncertain = rbeta (1,1,1)
      , lambda = rnorm(21,0,1)
      ,'.RNG.name'="base::Wichmann-Hill"
      ,'.RNG.seed'= this.seed
    )
  )
}

  date()
  set.seed (this.seed+1)

  compete.model <- run.jags(
   model=compete.bis.v,
   monitor=compete.parameters,
   n.chains=1,
   data=compete.data,
   inits=compete.inits.1(),
   thin=30, burnin=4000, sample=100,
   check.conv=FALSE, plots=FALSE)

  save.name <-  paste0 ("mexico_competing_final_result_sample_", which.sample,
                      which.process, ".Rda")

  chainsComp <- compete.model$mcmc[[1]]  
  # Uncomment next line to save files
  # save (chainsComp, file=save.name)
}
}

